#ifdef _OPENMP
#include <omp.h>
#endif
#include <iostream>
#include<stdlib.h>
#include <windows.h>
#include <nmmintrin.h>
#include <immintrin.h>
using namespace std;
const int N = 1000;
const int NUM_THREADS = 6;
float m[N][N];
int i, j, k;
float tmp;
void m_reset() {
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < i; j++)
			m[i][j] = 0;
		m[i][i] = 1.0;
		for (int j = i + 1; j < N; j++)
			m[i][j] = rand();
	}
	for (int k = 0; k < N; k++)
		for (int i = k + 1; i < N; i++)
			for (int j = 0; j < N; j++)
				m[i][j] += m[k][j];
}
void serial() {
	for (int k = 0; k < N; k++)
	{
		for (int j = k + 1; j < N; j++)
			m[k][j] = m[k][j] / m[k][k];
		m[k][k] = 1.0;
		for (int i = k + 1; i < N; i++)
		{
			for (int j = k + 1; j < N; j++)
				m[i][j] = m[i][j] - m[i][k] * m[k][j];
			m[i][k] = 0;
		}
	}

}
void omp_LU() {
#pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
	for (k = 0; k < N; ++k) {
		// 串行部分，也可以尝试并行化
#pragma omp single
		{
			tmp = m[k][k];
			for (j = k + 1; j < N; ++j) {
				m[k][j] = m[k][j] / tmp;
			}
			m[k][k] = 1.0;
		}
		// 并行部分，使用行划分
#pragma omp for
		for (i = k + 1; i < N; ++i) {
			tmp = m[i][k];
			for (j = k + 1; j < N; ++j) {
				m[i][j] = m[i][j] - tmp * m[k][j];
			}
			m[i][k] = 0.0;
		}
	}
}
void omp_SIMD() {
#pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
	for (k = 0; k < N; ++k) {
		int a = 0;
		// 串行部分，也可以尝试并行化
#pragma omp single
		{
			tmp = m[k][k];
			__m128 vt = _mm_loadu_ps(&m[k][k]);
			for (int j = k + 1; j + 4 < N; j += 4, a = j) {
				__m128 va = _mm_loadu_ps(&m[k][j]);
				va = _mm_div_ps(va, vt);
				_mm_store_ps(&m[k][j], va);
			}
			for (int j = a + 1; j < N; j++)
				m[k][j] = m[k][j] / m[k][k];
			m[k][k] = 1.0;
		}
		// 并行部分，使用行划分
#pragma omp for
		for (i = k + 1; i < N; ++i) {
			tmp = m[i][k];
			__m128 vaik = _mm_loadu_ps(&m[i][k]);
			for (int j = k + 1; j < N - 3; j += 4, a = j) {
				__m128 vakj = _mm_loadu_ps(&m[k][j]);
				__m128 vaij = _mm_loadu_ps(&m[i][j]);
				__m128 vx = _mm_mul_ps(vakj, vaik);
				vaij = _mm_sub_ps(vaij, vx);
				_mm_store_ps(&m[i][j], vaij);
			}
			for (int j = a + 1; j < N; j++) {
				m[i][j] = m[i][j] - m[k][j] * tmp;
			}
			m[i][k] = 0.0;
		}
	}
}
int main() {

	m_reset();
	long long head, tail, freq;
	QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	serial();
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "LU分解单线程:" << (tail - head) * 1000.0 / freq << "ms" << endl;
	// 在外循环之外创建线程，避免线程反复创建销毁，注意共享变量和私有变量的设置

	m_reset();
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	omp_LU();
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "LU分解多线程:" << (tail - head) * 1000.0 / freq << "ms" << endl;
	

	m_reset();
	QueryPerformanceCounter((LARGE_INTEGER*)&head);
	omp_SIMD();
	QueryPerformanceCounter((LARGE_INTEGER*)&tail);
	cout << "SIMD多线程:" << (tail - head) * 1000.0 / freq << "ms" << endl;
//for (int i = 0; i < N; i++) {
//		for (int j = 0; j < N; j++)
//			cout << m[i][j]<<" ";
//		cout << endl;
//	}

}